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Abstract. Ornstein-Uhlenbeck processes driven by general Levy process are 
considered in this paper. We derive strongly consistent estimators for the 
moments of the underlying Levy process and for the mean reverting parameter 
of a discretely observed Levy driven Ornstein-Uhlenbeck process. Moreover, 
we prove that the estimators are asymptotically normal. We use ergodicity 
(*^ , arguments. Finally, we test the empirical performance of our estimators in a 

■ 4-^ ■ simulation study and we fit the model to real VIX data. 
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1. Introduction 

i^ ' Given a positive number A and a time- homogeneous Levy process L, the Ornstein- 

fvj . Uhlenbeck (OU) process driven by L is defined by 



> 

(N 



O ■ (1) Yt = e-^'Yo + e-^* / e'dL, 

OO ■ 

^D • where Yq is assumed to be independent of {Lt}t>o- Following the terminology 

introduced by Barndorff-Nielsen and Shephard in [2] , we shall call L the background 
driving Levy process (BDLP). It is easy to see that ((1} is the unique strong solution 
of the stochastic differential equation 






(2) dVt = -XYtdt + dLxf 

Under some regularity conditions on the Levy measure of L and if A > 0, y admits 
a unique invariant distribution Fy ■ Owing to the scaling of the time index of L in 
([2]) by A (i.e. the term L\t), Fy is independent of A. 

Let us suppose now that we have discrete-time observations Yo,Yh, ■ ■ ■ , Y(„_i)^ 
with h > from {Yt}t>o as it is defined by ([1]). The objective here is to estimate 
the parameters of the model using these discrete-time observations. In particular, 
we are interested in estimating A and moments of Li. We derive strongly consistent 
method of moments estimators and prove that they are asymptotically normal. In 
this paper, we consider A and only the first two moments of Li, i.e. ELi and ELf . 
However, the methodology can be extended to higher moments as well (see Remark 
13. 2p . Similar methods of moments estimators have been used elsewhere as well (see 
Valdivieso et al ^25, for example), but without theoretical justification. The proof 
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of their consistency and asymptotic normality is presented, according to the best 
knowledge of the author, for the first time in the present work. 

Our motivation for studying this problem comes from continuous stochastic 
volatility models in financial mathematics. BarndorfF-Nielsen and Shephard (in 
[2]; see also [3]) model stock price as a geometric Brownian motion and the diffu- 
sion coefficient of this motion as an OU process that is driven by a subordinator (a 
Levy process that is nonnegative and nondecreasing) . Other continuous stochastic 
volatility models can be found in Kliippelberg et al [15, and in Shephard [23]. Some 
papers that consider statistical inference of these models are Barndorff-Nielsen and 
Shephard \2., Brockwell et al i5_, Haug et al [12], Jongbloed et al [13] . 

The paper is organized as follows. Section 2 presents several known results 
concerning OU processes and Levy processes. In section 3 we consider strongly 
consistent estimators of the first two moments of Li and of A and we provide a 
methodology to express any moment of the stationary distribution of {Y}t>Q in 
terms of the moments of Li. In section 4, we prove that these estimators are 
asymptotically normal. Section 5 discusses modeling issues and simulation tech- 
niques and presents simulation results for gamma OU process and inverse Gaussian 
OU process. In section 6 we fit the model to real log(VIX) data and we argue that 
an OU model is a good candidate for modeling log(VIX). Finally, section 7 contains 
a summary and a discussion on future work. 

We would like to mention here, that after completion of this work, the author 
learned about the results in Jongbloed et al 13 . In [13 , the authors assume that L 
is a subordinator. Let _FY denote the Levy measure of L, Y the unique stationary 
solution to ([Ij (which exists if J^ ^\og{x)F]^{dx) < oo for example) and Fy its 
probability law. The characteristic function of Y is given by 

^p^(t) := j e'^'^Fyidx) = exp( /°°[e**" - l]'^dx), 



J JQ X 

where k{x) — Fl{x, cx)). Hence, the stationary distribution Fy , of the OU process 
Y, is being determined by the canonical function Hi(x). In [13j . the authors develop 
a nonparametric inference procedure for A and for the canonical function k{x). The 
results in the present paper complement the results of |13] . 

2. Assumptions and Preliminary Results 
Consider a probability space (51,5,^) equipped with a filtration 5^^. 

Definition 2.1. A one dimensional ^t adapted Levy process is usually denoted by 
Lt — Lt{uj), t > 0, u! G fl and is a stochastic process that satisfies the following: 

(i) Lt e dt for all t > 0. 

(ii) Lq ^ a.s. 

(iii) Lt — Ls is independent of ^s o,nd has the same distribution as Lt-s- 
(iv) It is a process continuous in probability. 

We assume that we are working with a cadlag Levy process (i.e. it is right 
continuous with left limits). It is well known that every Levy process has such a 
modification. 

Furthermore, if F^ denotes the Levy measure of Li, we will assume that there 
exist a constant A/ > such that 

(3) ( e'^^FLidx) < oo, for every \v\ < M. 

J\x\>l 
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Condition ([3]) gaurantees that the moment generating function v — > Ee"'^^ exists at 
least for |t)| < M (see Wolfe ^ and Eberlein and Raible [9]). 
We shall write 
(i) ELi = ^. 
(ii) Var(Li) ^ a^. 
Moreover, we shall assume that Yq is independent of {Lt}t>o and that 

(4) Yo ^ 

/o 

The integral on the right hand side of (|3]) is well defined (see Sato [35] for example). 

The following proposition, which is a reformulation of Propositions 1 and 2 in 

Brockwell 6 , characterizes the stationarity of the OU process {Yt}t>o- 

Proposition 2.2. If Yq is independent of {Lt}t>o and EL\ < cx) then {Yt}t>o is 
weakly stationary if and only if X > and Yq has the same mean and variance as 
/p e~^dLs. If in addition Yq has the same distribution as /g e~^dLs, then {Yt}t>o 
is strictly stationary and vice-versa. 

In Masuda [I9j now, the author proves, under mild regularity conditions, that the 
OU process Y is strong Feller, its probability law has a smooth transition density, 
is ergodic and exponentially /3— mixing (strong mixing). Before mentioning the 
results of [TO] that we will use in the present paper, let us recall the definitions of 
a self-decomposable law on R and of /?— mixing. 

Definition 2.3. Let X be a positive number. Then, an infinitely divisible distribu- 
tion Fy is called X— self-decomposable, if there exists a random variable X — Xt^\, 
such that, for each t G M+ 

4>Fy{u) — </'Fy(e~'^*u)0Fx ("))"" G ^: 

where (Jjfy ("") ^'^'^ 4'Fx (jJ-) '^''g i^^ characteristic functions corresponding to Fy and 
Fx respectively. For the sake of notational convenience we will just say that Fy is 
called self-decomposable. 

If Ji^i^i \og{\x\) FL{dx) < cx), then the class of all possible invariant distributions 
of Y forms the class of all self-decomposable distributions Fy (see Sato [H]). In 
particular, the latter is implied by ([3]). 

Definition 2.4. For a stationary process Y = {Yt}t>o define the a-algebras ^i — 
S'CCu) = o-{{Yy},0 <v<u) and ^2 = ^{u+t,oo) = o-{{Yv},v >u + t). Then 
(i) Y is called j3— mixing (or strong mixing) if: 

j3{t)^ sup \P{An B) - P{A)P{B)\ ^0 ast^ oo. 

(ii) Y is called /3— mixing with exponential rate if for some fc > and a > 0: 

(3{t) < ke~°-* for t > 0. 

The following theorem is Theorem 4.3 in Masuda [19] and discusses the mixing 
properties of {Yt}t>Q. 

Theorem 2.5. Let A > and {Yt}t>o be the strictly stationary OU process given 
by (QP with self-decomposable marginal distribution Fy . If we have that 

\xfFydx < oo 
3 
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for some p > 0, then there exists a constant a > such that /3{t) ~ 0(e^°*) as 
i — ?> cx). In particular, Y is ergodic. 

3. Method of Moments Estimation 

We aim at estimation of the model parameters Oq — (/i,cr^, A) from a sample of 
equally spaced observations from ([T]) by matching moments and empirical autocor- 
relation function to their theoretical counterparts. 

Proposition 13. II below relates the theoretical moments of ii with the theoretical 
moments of the stationary distribution Fy of {Yt}. 

Proposition 3.1. Suppose that {Lt}t>o is a Levy process such that ELi — fi < oo, 
VarLi = a^ < oo and that (0) holds. Let M be the largest constant satisfying (0) 
and assume that X < M . Then, the following are true 

(i) EFo = A* 
(ii) VarYo = 2l 

Proof. Let j(v) be the cumulant function of Li, i.e. 

(5) 7(w) =lnEe"-^i 

By the Levy- Khinchine representation Theorem we get that "f{v) has the form 

(6) 7(i;) = bv+^v^+ f (e^^ - 1 - vx)FL{dx), 

which is valid for \v\ < M . Moreover, 7 is continuously differentiable (see Lukacs 

m)- 

Using the assumptions ELi = /i and VarLi ~ a^ and relations ([5]) and ([5]), it is 
easy to see that h — ^i and c = ct^ — /^ x^F^idx). 

In order to calculate EYo and VarYo we use the following formula: 

(7) Ee-'o"^''"''^-^-' ^(^ir'ri^^^")''-^^ 

which is valid since X < AI (see Lemma 3.1 of Eberlein and Raible j^). 

Recall now that we have assumed Yq = J^ e^^dLg in distribution. The latter 
and (O imply that: 











EYo 


= -^Eelo 
dX 

d r^ 

dX 


7(Ae-= 


|a=o 

)ds| 


(8) 










= M 






In 


a similar 


way 


we 


get that 


EY,^ = 4 


+ M^ 


. This 


pre 


iposition. 















This concludes the proof of the 

D 



Remark 3.2. We would like to note here, that the proof of Proposition \3A\ can be 
used for the calculation of higher moments of Yq . 

It follows directly by ([T]) that the theoretical autocovariance and autocorrelation 
function of Yt are given by the formulas 

(i) autocovariance: 7(/i) = cov{Yt+h,Yt) — ^e^^'', for h e Nq. 
(ii) autocorrelation: p{h) = corr{Yt+h,yt) — e^^'^, for h G Nq. 
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On the other hand, the empirical moments, autocorrelation and autocovariance 
function are given by the formulas below. Let d > be fixed. Then, we have: 

(i) Sample mean: Y. = ^ J27=i ^i- _ 
(ii) Sample variance: ^ X]r=i(^' ^ ^)^- 
(iii) Sample autocovariance: 7„ = (7„(0),7„(1), • • • ,7„(d))"^ where for h G 

{0, • • • , 4 we define %ih) = i E7=iiY,+h - Y.)iY, - Y). 
(iv) Sample autocorrelation: /5„ = (/5„(0), p„(l), • • • , pn{d))'^ where for h E 
{0, • • • ,d} we define p„(/i) = ?^. 
We have the following Theorem: 

Theorem 3.3. Let ^, (t^,7(-),7'„(-), p(-) and pn{-) be defined as above. Then, the 
following statements are true 

(i) Y. — > ji almost surely 

(ii) i Er=i(^^ - y-? "-^ ^ <^^rnost surely 

(iii) (7„(1), • • • ,7n(d)) "-^ (7(1)^ • ■ ■ ^l{d)) almost surely 

(iv) (/5„(1), • ■ ■ , Pn{d)) "-^ (p(l), • • • , p{d)) almost surely 

Proof. Due to our assumptions, the process {Yt\t>a is strictly stationary. Moreover 
by Theorem 12.51 it is also /3— mixing with exponential decaying rate. These two re- 
sults imply ergodicity of {Yt}t>o. The latter together with strict stationarity imply 
that empirical moments and sample autocovariance functions are strongly consis- 
tent estimators of the corresponding theoretical quantities Billingsley [4 . Then, 
the statement of the Theorem follows. D 

For the mean reverting parameter A we have the following Lemma. 

Lemma 3.4. Let K be a compact subset of M+ such that the true value of \, say 
Xo, belongs to K and let Xn = fl1?"^*'^AG_fs: S/i=i(Pn(^) ^ e^'*''')^. Then A„ exists, 
is locally unique and 

(9) A„ — > Xo almost surely. 



Proof Consider the functions A„ (A) = Y.h=iiMh)-P\[h)Y andAo(A) = I]/j=i(pa„(/i)- 
p\{h)Y where p\{h) = e~^'^ . Theorem [331 implies that for all Xe K: 

A„(A) — > Ao(A) almost surely. 

By Theorem ILl in Andresen and Grill [T we have 

sup |A„(A) — Ao(A)| — > almost surely. 

Observe now that Ao(A) is a sum of nonnegative terms. It becomes zero if and only 
if A = Aq. Hence, Ao(A) has a unique minimum at A = Aq which is equal to zero. 
We get 

A„(Ao) — > almost surely. 
Furthermore, for n finite we have that < A„(A„) < A„(Ao). Therefore we get 

A„(A„) — > almost surely. 
5 
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Moreover, we have 

d 

|A„(A„)-Ao(A„)| = \Y,{plih)-pUh) + 2p^^^^{h){p^,{h)~p,M)]\< 

h=l 
d 

< T.^\Pn{h)\ + IpAo Wl + 2\p^Jh)\]\p^,{h) - p„(/i)| < 
h=l 

d 

< 4yj \p\g{h) — Pn(h)\ ""- — "^ almost surely. 

Here, we used the relation |p„(ft,)| < 1 which follows immediately from Cauchy- 
Schwarz inequality. The above imply that 

Ao(A„) — > almost surely. 

But, Aq is the unique minimum of Ao(A) and it satisfies Ao(Ao) — 0. Thus, 

Ao(A„) " — ^ Ao(Ao) — almost surely. 

Hence, we easily conclude (Corollary H.2 in [1]) that An is locally uniquely deter- 
mined and that 

A„ — > Xo almost surely. 

D 



Remark 3.5. Theorem \3.3\ and Lemma[3^ give us two strongly consistent esti- 
mators for A. The first one is Ai,n — — log(/5„(l)) and the second one is A2,n = 
argmin)^^f^^-^{pn(h) ~ e^'^^Y' . One could use, for example, Ai,„ as an initial value 
to an algorithm that calculates \2,n- 

Summarizing, we have that /i„, tr^ and Ai,„, \2,n are strongly consistent estima- 
tors of /i, (T^ and A respectively, where: 



(10) 



An = 




-I - 


- 2l>;(r,-An)' 

n ^-^ 


Al,n = 


= -log(p„(l)) 


A2,n = 


= argmin;,)^(p„(/i)-e ^'') 



h=l 

Remark 3.6. For a stationary model, the parameter A has to be positive. However, 
if we compute X2.n CiS the unrestricted minimum, A2,n = argmin-y^^ S/i=i(Pn(^) ~ 
g-A/i-j2 ^g rnay end up with a negative estimator A„. In this case, we define the 
estimator of A to be zero and we take this as an indication that the data is not 
stationary. 
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4. Asymptotic Properties of the Moment Estimators 



In this section we prove that the estimators defined by (jlO[) are asymptotically 
normal. 

If /? is a vector, then we define by /3'^ its transpose. We begin with the following 
central limit theorem. 

Theorem 4.1. Let us assume that there exists a S > such that EYq'^ < oo. 
Define 

tpn = (/i„,7„(0),7„(l),--- ,7„(d))^ 
V'o = (Ai,7(0),7(l),---,7(rf)r 
S — [<^fc./]fcl^i with elements 

oo 

aku = cov{Z'^, Z{) + 2^ cov{Z'^, Z-_^ J where 

i=l 

Then, the following holds: 

(11) V^(V3„ - ^o) A iV(0, E) 

where N{0, S) is t/ie multivariate normal distribution with mean anrf variance- 
covariance matrix S. 

Proof. The proof of this theorem is similar to the proof of Proposition 3.7 of Haug 
et al [12]. Let us define 

(i) iiih) = ^Y:^=l(y^+h - M)m - m),/^ e {0, • • • , di- 
al) 7: = (7:(o),---,7:(rf))^. 

We first prove that ([11]) is true with f/)* = (/i„,7*(0),7*(l), • • • ,7^(d))^ in place 
of -0„. 

By the well known Cramer- Wold device, it is sufficient to prove that for every 
/3 e K'^+2 such that /3'^E/3 > we have 



1 " 
(12) V^(-V/3^Z,-/3^Vo) AiV(0,/3^S/?). 

Tt ^ ^ 



n 



It is well known (see [U for example) that strong mixing and the corresponding de- 
caying rate are preserved under linear transformations. Thus, the sequence {f3'^ Zi} 
is strong mixing with exponential decaying rate. Since, by assumption E|Zip+'^ 
for some e > 0, the central limit theorem for strong mixing processes is applicable 
(Theorem 7.3.1 in Ethier and Kurtz TO] ) . Hence, we have as n — > oo that 

n 

(13) v^(-V/3^Z,-/3^Vo) AiV(0,a2). 

n ^-^ 

i=\ 

But, we easily see that ct^ = var{0^ Z{) + 2X;^i cou(/3^ Zi,/3'^Zj+i) = /3'^E/3. So 
dm holds. 

Now recall that by Theorem 13.31 we have 



(14) ■!/;„ — > Vo almost surely. 

7 
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Following the proof of proposition 7.3.4 of Brockwell and Davis [^ '^^ S^t 

1 " 
(15) y^(- y /3^Z, - /3^V>„) "-^ in probability. 

1—1 

Therefore, V'n has the same asymptotic behavior as V'^. The latter and (J12p imply 
the Theorem. D 

Corollary 4.2. Let the conditions of Theorem \4.1\ hold. Then we have 

^^(p„-p) AiV(0,I]p). 

Proof. It follows directly by Theorem l4.1l and delta method (Theorem 3.1 in A.W.van 
der Vaart [IS]). D 

Finally, we prove central limit theorem for 0„ = {fin,(j'^, X2,n)'^ ■ Let us denote 
CTy = 7(0) = Var(yo) = %- and define the following mappings. 



(16) G : M X [0, 00)^ — >Rx [0, ooY : G{n, af , A) = 



(Ai,2crf,,A), A>0 
(Ai,24,0), A<0. 



d 



(17) F : R^+i -^ M+ : ^(/}) = argmin;, ^(/5«(/i) - e-^'')^ = A2,„ 

and iJ as follows: 

(18) H : M''+2 — ^ R X [0, oo)^ : H{p, 7^) = G(^, crf^, F(p)), 
where p{h) = ^4^1 for /i = 0, • • • , d. 

Theorem 4.3. Lei t/ie conditions of Theorem \4. 1\ be satisfied. Let us define: 

^n = (Aril O'n, A2,„) 

Then the following holds: 

(19) ^^(^„-0o)^[^i^^]iV(O,E) 

9(^,7^) 

Proof. It follows directly by Theorem 14.11 and delta method applied to the differ- 
entiable map H. D 

5. Modeling and Simulation 

In this section we discuss modeling issues of Levy driven OU processes and 
present some simulation results for a gamma OU process and an inverse Gaussian 
OU process. We use the simulated data to test the performance of our estimators. 
However, first we mention the necessary ingredients for the simulation process. 

A very important ingredient in modeling of Levy driven OU processes is the 
connection between the Levy density of the stationary distribution of Y to the 
Levy density of the probability law of Li. In particular we have the following 
proposition. 
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Proposition 5.1. Assume that the Levy density oJY, uy{x), is differ entiahle and 
denote the Levy density of the probability law of Li by V]^{x). Then the following 
relation holds. 

(20) fhix) = —vy{x) — xv'y{x). 



Proof. It follows directly by the fact that the stationary solution, F, to ^ satisfies 

/•oo 

yi / e-^'dL{\s). 







See P] and [3] for more details. D 

Hence, given Vhix) we can find vy{x) and vice- versa. One can specify the law of 
the one dimensional marginal distribution of the OU process Y and work out the 
density of the BDLP, Li. One can also go the other way and model through the 
BDLP. Of course, there are constraints on valid BDLP's which must be satisfied. 
In particular, if 

min{l,a; ^VL{x)dx 



then vl{x) is the density of a Levy jump process L and there exists an OU process 
Y such that L is the BDLP of Y . A very good survey on the relation between 
several distributions of Y and L is BarndorfF-Nielsen and Shephard JS] (see also 
Schoutens [25]). 

Another important ingredient in simulations is the infinite series representation 
of Levy integrals (see Rosinski ^21;). For simplicity, we restrict attention to Levy 
processes, L, that are subordinators, i.e. they are nonnegative and nondecreasing. 
It is easy to see that subordinators have no Gaussian component, nonnegative drift 
and a Levy measure that is zero on the negative half-line. If Y models stochastic 
volatility then it has to be positive and such a choice of the BDLP guarantees that. 

Let us denote by Pj the tail mass function of v^, i.e. 



rj(x) = / 

J X 



(21) r+(x) = / VL{y)dy 

J X 

and by F^ the generalized inverse function of Pj, i.e. 

(22) Tl\x)^ini{y>Q:Tl{y)<x}. 



rXt 



In order to simulate from ([T]) we need to be able to simulate from e~'^* L e'^dLg 
The key result here is the following infinite series representation of this type of 
integrals (Rosinski |21j): 

Proposition 5.2. Consider a subordinator L with positive increments. Let f be a 
positive and integrable function on [0,r]. Then 

(23) / f{s)dLs^Y.^l\a,lT)f{Tri), 

where the equality is understood in distributional sense, {cm} and {ri\ are two 
independent sequences of random variables such that ri are independent copies of 
a uniform random variable in [0, 1] and {oti\ is a strictly increasing sequence of 
arrival times of a Poisson process with intensity 1 . 

9 
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Remark 5.3. 

slow. 



We note here that the convergence of the series 



is often quite 



Using (|23|) we can then simulate a Levy driven OU process. In particular, if A 
denotes the time step, we will use the identity 

rt+A 



Yt 



t+A 



= e-^^(r, 



,-At 



'dLxs) 



(24) 



-AA 



(Yt 



'dLxs). 



Let us demonstrate the validity of our estimators modeling through the BDLP. We 
consider two cases: (a) when Yq ^ Gamma(a, b) and (b) when Yq '^ IG(a, b), where 
IG stands for inverse Gaussian. 

Regarding the A parameter, we recall our estimators: Ai^n = — -2S1^-LJ1 and 

A2,n = argmiu;^ J2h=i(Pn{h) — e"'^'"^)^. In ([T0|) we defined Ai,„ and A2,n for A = 1, 
but of course one can generalize them to any A > 0. 

5.1. Gamma OU modeL Assume that the driving Levy process L is a compound 
Poisson process and in particular, that Lt = J2nLi ^n where Nt is Poisson with 
intensity parameter a and Xn are independent identically distributed Gamma(l,6) 
random variables. Using (|20l) we get that Yq ~ Gamma(a, 6). It is known (see [2]) 
that in this case 

r£i(x)=max{0,-ilog(-)}. 
a 

Using this and equations ([^^ and (IMl) we can easily simulate from a Gamma(a, 6)- 
OU process. We also need to know how the parameters fj, and a^ relate to a and b. 
Since ELi = fi and Var(Li) — a^ implies EYq — ^ and Var(lo) = ^, we have that 

a = 24 and 6 = 24- 

We simulated 100 independent paths of a gamma OU process of 1000 observa- 
tions each, with time step A = 0.1, using (p4)) . We chose fj, = 2, a^ = 0.25 and 
in order to capture possible different behaviors of the intensity parameter we chose 
two different values for A, 0.5 and 5. 

Tables I and II, summarize the results for Oq = (2, 0.25, 0.5) and for Oq = 
(2, 0.25, 5) respectively. 



True Values 


Est. Values 


Sample Std. Error 


Comments 


M = 2 


1.995458 


0.0702198 


- 


cr^ = 0.25 


0.2350207 


0.05352894 


- 


A = 0.5 


0.566116 


0.1126439 


An = Ai^n 


A = 0.5 


0.5879571 


0.1441501 


An = A2,„ 



Table 1. 6*0 = (2,0.25,0.5) 
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True Values 


Est. Values 


Sample Std. Error 


Comments 


M = 2 


2.003799 


0.02094129 


- 


a' = 0.25 


0.2473567 


0.01608991 


- 


A = 5 


5.12962 


0.4463517 


An = Ai^n 


A = 5 


5.186585 


0.5898125 


An = ^2,n 



Table 2. 6*0 = (2,0.25,5) 



5.2. Inverse Gaussian OU model. It is well known that if Yq 
the Levy density of Y is 

(25) i^y{x) 



IG{a,b), then 



/2n 



-3/2^~^b^x_ 



Consider now the Lambert-W function, Lwi-), which satisfies Ltu(a;)e^™('^) = x. As 
it is also shown in Gander and Stephens [IT], equations (PH)) . (pij) and (1^ imply 
that the inverse tail mass function of the BDLP of an IG(a,b)-OU process is given 

by 



(26) 



1 



a%^ 



r-^Hx) - T^L^ii^-:^) 



62 '"^27rx2' 

Using the latter and equations ([25]) and ([M|) we can easily simulate an IG(a,b)-OU 
process. We also need to know how the parameters /z and a'^ relate to a and b. 
Since ELi = fi and Var{Li) = a" implies EYq — /i and Var{Yo) = ^, we have 

that a — u\/ ^ and b — ■ ' — 



We simulated 100 independent paths of an IG-OU process of 1000 observations 
each with time step A = 0.1, using ^^. As before, we chose fi — 2, cr^ ~ 0.25 and 
two different values for A, 0.5 and 5. 

Tables III and IV, summarize the results for 6*0 = (2, 0.25, 0.5) and for 6*0 = 
(2, 0.25, 5) respectively. 



True Values 


Est. Values 


Sample Std. Error 


Comments 


/i = 2 


1.986862 


0.06476202 


- 


0-2 = 0.25 


0.2331244 


0.05235387 


- 


A = 0.5 


0.5581237 


0.1128397 


An = Ai.n 


A = 0.5 


0.6050457 


0.1376689 


An = A2.„ 



Table 3. 0o = (2,0.25,0.5) 



Table 4. 



(2,0.25,5) 



True Values 


Est. Values 


Sample Std. Error 


Comments 


/. = 2 


1.955288 


0.03107831 


- 


a' = 0.25 


0.2452349 


0.01750871 


- 


A = 5 


5.05211 


0.4262788 


An = M,n 


A = 5 


5.158421 


0.659508 


An = A2,n 
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6. Real Data Analysis 

In 1993, the Chicago Board Options Exchange (CBOE) introduced the CBOE 
volatihty index, VIX, and it quickly became a popular measure for stock market 
volatility. In 2003, the VIX methodology was updated (see www.cboe.com for 
more details on the old and new VIX methodology). VIX measures the implied 
volatility of S&P 500 index options and it provides a minute-by-minute snapshot 
of the markets expectancy of volatility over the next 30 calendar days. 

We fitted the gamma OU model and the IG OU model to daily log opening 
values of the VIX for the year 2004 (VIX values are calculated using the new 
methodology). The data are taken from www.cboe.com. We use the values from 
1/2/2004 tin 9/30/2004 for the calibration of the model (in total 189 data points) 
and the values from 10/1/2004 tiU 11/30/2004 (in total 41 data points) for testing 
the model. 

Table V summarizes the estimators, given by (jlOp . of the parameters of the 
model. We used A2 n to estimate A. 



parameter 


M 


a' 


A 


estimated value 


2.781769 


0.01919740 


0.1767250 



Table 5. Estimated values for the parameters. 



In Figure 1, we see the first 10 lags of the empirical autocorrelation function 
of the log(VIX) for 1/2/2004 tiU 9/30/2004 versus the theoretical autocorrelation 
function of the OU model with A = 0.1767250, i.e. p{h) = e-oi767250^^ 















rrr ™5j;™pT°3°ror' 

































Figure 1 . True time series acf versus the model implied acf with 
A = 0.1767250. 



As we saw before, the autocorrelation function of an OU model is exponentially 
decreasing, i.e. it has the form e^'^'*. Figure 1 shows that e-oi767250h approximates 
sufficiently well the empirical autocorrelation function of log(VIX) for 1/2/2004 till 
9/30/2004, which is also exponentially decreasing. Hence, we conclude that an OU 
model is a good candidate for describing this data set. 

To investigate the model fit, we performed a Ljung-Box test for the squared 
residuals. We used the estimated values from Table V and since our data is daily 
opening values we chose A = 1 for the time step. The test statistic used 10 lags 
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of the empirical autocorrelation function. The gamma OU model performed better 
than the the IG-OU model. The null hypothesis was not rejected at the 0.05 
level and the p— value was quite high, 0.6165. In Figure 2 we see the empirical 
autocorrelation function of the residuals of log(VIX) and in Figure 3 we see the 
actual residuals of the gamma OU model. 

In Figure 4 we see in one figure: the actual time series from 10/1/2004 till 
11/30/2004, the one step ahead predicted time series and 95% bootstrap upper 
and lower confidence bounds of the one step ahead predicted time series. In order 
to create the one step ahead predicted time series we averaged over 50 paths. We 
observe that the real time series (solid line) is most of the time within the 95% 
bootstrap upper and lower confidence bounds of the one step ahead predicted time 
series (dotted lines), with very few exceptions. 



Figure 2. Empirical acf for the residuals of log(VIX). 



Figure 3. The actual residuals. 
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Figure 4. Actual time series versus one step ahead predicted time series. 



7. Discussion And Future Work 

In this paper, we consider an Ornstein-Uhlenbeck process driven by a general 
Levy process. We derive strong consistent estimators for the parameters of the 
model and we prove that they are asymptotically normal. Using simulated data, 
we show that the estimators perform well at least for a gamma OU model and for 
an IG-OU model. Lastly we fit the model to real data and we see that a Levy 
driven OU model is a good candidate for describing log(VIX). 

There are some interesting extensions to the model studied in this paper. One 
such extension is a coupled two dimensional OU process driven by a two dimensional 
Levy process. This model is important for financial applications, since it could 
be used to model log of the price and stochastic volatility simultaneously. An 
interesting question for financial applications is option pricing in these type of 
models (see Nicolato and Venardos [20 for some recent related results). These 
questions will be addressed in future work. 
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